{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0b149598",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "OMP: Info #270: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1  3  9 27 30 35]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import galois\n",
    "\n",
    "# R1CS to QAP\n",
    "\n",
    "GF = galois.GF(67)\n",
    "\n",
    "A = np.array([\n",
    "    [0, 1, 0, 0, 0, 0],  # 约束 1\n",
    "    [0, 0, 1, 0, 0, 0],  # 约束 2\n",
    "    [0, 1, 0, 1, 0, 0],  # 约束 3\n",
    "    [5, 0, 0, 0, 1, 0],  # 约束 4\n",
    "])\n",
    "\n",
    "\n",
    "B = np.array([\n",
    "    [0, 1, 0, 0, 0, 0],  # 约束 1\n",
    "    [0, 1, 0, 0, 0, 0],  # 约束 2\n",
    "    [1, 0, 0, 0, 0, 0],  # 约束 3\n",
    "    [1, 0, 0, 0, 0, 0],  # 约束 4\n",
    "])\n",
    "\n",
    "C = np.array([\n",
    "    [0, 0, 1, 0, 0, 0],  # 约束 1\n",
    "    [0, 0, 0, 1, 0, 0],  # 约束 2\n",
    "    [0, 0, 0, 0, 1, 0],  # 约束 3\n",
    "    [0, 0, 0, 0, 0, 1],  # 约束 4\n",
    "])\n",
    "\n",
    "A_galois = GF(A)\n",
    "B_galois = GF(B)\n",
    "C_galois = GF(C)\n",
    "\n",
    "x = GF(3)\n",
    "v1 = x * x\n",
    "v2 = v1 * x\n",
    "v3 = v2 + x\n",
    "y = v3 + GF(5)    \n",
    "\n",
    "witness = GF(np.array([1, x, v1, v2, v3, y]))\n",
    "\n",
    "print(witness)\n",
    "\n",
    "assert all(np.equal(np.matmul(A_galois, witness) * np.matmul(B_galois, witness), np.matmul(C_galois, witness))), \"not equal\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "08c84494",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "poly A:  [Poly(12x^3 + 62x^2 + 65x + 62, GF(67))\n",
      " Poly(44x^3 + 5x^2 + 11x + 8, GF(67))\n",
      " Poly(34x^3 + 63x^2 + 43x + 61, GF(67))\n",
      " Poly(33x^3 + 37x^2 + 60x + 4, GF(67))\n",
      " Poly(56x^3 + 66x^2 + 13x + 66, GF(67)) Poly(0, GF(67))]\n",
      "poly B:  [Poly(22x^3 + 36x^2 + 6x + 3, GF(67))\n",
      " Poly(45x^3 + 31x^2 + 61x + 65, GF(67)) Poly(0, GF(67)) Poly(0, GF(67))\n",
      " Poly(0, GF(67)) Poly(0, GF(67))]\n",
      "poly C:  [Poly(0, GF(67)) Poly(0, GF(67)) Poly(11x^3 + 35x^2 + 18x + 4, GF(67))\n",
      " Poly(34x^3 + 63x^2 + 43x + 61, GF(67))\n",
      " Poly(33x^3 + 37x^2 + 60x + 4, GF(67))\n",
      " Poly(56x^3 + 66x^2 + 13x + 66, GF(67))]\n"
     ]
    }
   ],
   "source": [
    "## low-degree extension for A, B, C\n",
    "def interpolate_column(col):\n",
    "    xs = GF(np.array([1,2,3,4]))\n",
    "    return galois.lagrange_poly(xs, col)\n",
    "\n",
    "# axis 0 is the columns. apply_along_axis is the same as doing a for loop over the columns and collecting the results in an array\n",
    "A_polys = np.apply_along_axis(interpolate_column, 0, A_galois)\n",
    "B_polys = np.apply_along_axis(interpolate_column, 0, B_galois)\n",
    "C_polys = np.apply_along_axis(interpolate_column, 0, C_galois)\n",
    "\n",
    "# matrix A becomes 6 polynomials\n",
    "print(\"poly A: \", A_polys)\n",
    "print(\"poly B: \", B_polys)\n",
    "print(\"poly C: \", C_polys)\n",
    "\n",
    "# poly A:  \n",
    "#  [Poly(12x^3 + 62x^2 + 65x + 62, GF(67))\n",
    "#  Poly(44x^3 + 5x^2 + 11x + 8, GF(67))\n",
    "#  Poly(34x^3 + 63x^2 + 43x + 61, GF(67))\n",
    "#  Poly(33x^3 + 37x^2 + 60x + 4, GF(67))\n",
    "#  Poly(56x^3 + 66x^2 + 13x + 66, GF(67)) \n",
    "#  Poly(0, GF(67))]\n",
    "# "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e1918ac5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4x^6 + 18x^5 + 3x^4 + 13x^3 + 38x^2 + 12x + 46\n"
     ]
    }
   ],
   "source": [
    "## get A B C term\n",
    "from functools import reduce\n",
    "\n",
    "def inner_product_polynomials_with_witness(polys, witness):\n",
    "    mul_ = lambda x, y: x * y\n",
    "    sum_ = lambda x, y: x + y\n",
    "    return reduce(sum_, map(mul_, polys, witness))\n",
    "\n",
    "A_term = inner_product_polynomials_with_witness(A_polys, witness)\n",
    "B_term = inner_product_polynomials_with_witness(B_polys, witness)\n",
    "C_term = inner_product_polynomials_with_witness(C_polys, witness)\n",
    "\n",
    "print(A_term * B_term - C_term)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "03ff720d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4x^6 + 18x^5 + 3x^4 + 13x^3 + 38x^2 + 12x + 46\n"
     ]
    }
   ],
   "source": [
    "## Target polynomial T = (x-1)(x-2)(x-3)(x-4)\n",
    "T = galois.Poly([1, 67-1], field = GF) * galois.Poly([1, 67-2], field = GF) * galois.Poly([1, 67-3], field = GF) * galois.Poly([1, 67-4], field = GF)\n",
    "\n",
    "## Quotient polynomial Q = (A * B - C) / T\n",
    "Q = (A_term * B_term - C_term) // T\n",
    "\n",
    "print(T * Q)\n",
    "\n",
    "poly_big = A_term * B_term - C_term\n",
    "assert poly_big // T * T == poly_big"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f013266a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4x^6 + 18x^5 + 3x^4 + 24x^3 + 39x^2 + 66x + 47\n"
     ]
    },
    {
     "ename": "AssertionError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m/var/folders/xz/sj6gs0f150n67pbx3mg32v3c0000gn/T/ipykernel_65295/4053868161.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0mpoly_big_wrong\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mA_term_wrong\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mB_term_wrong\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mC_term_wrong\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m \u001b[0;32massert\u001b[0m \u001b[0mpoly_big_wrong\u001b[0m \u001b[0;34m//\u001b[0m \u001b[0mT\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mT\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mpoly_big_wrong\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     13\u001b[0m \u001b[0;31m## Error: A*B - C can not divide T\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAssertionError\u001b[0m: "
     ]
    }
   ],
   "source": [
    "## wrong witness\n",
    "\n",
    "witness_wrong = GF(np.array([1, 3, 9, 27, 30, 36]))\n",
    "\n",
    "A_term_wrong = inner_product_polynomials_with_witness(A_polys, witness_wrong)\n",
    "B_term_wrong = inner_product_polynomials_with_witness(B_polys, witness_wrong)\n",
    "C_term_wrong = inner_product_polynomials_with_witness(C_polys, witness_wrong)\n",
    "\n",
    "print(A_term_wrong * B_term_wrong - C_term_wrong)\n",
    "\n",
    "poly_big_wrong = A_term_wrong * B_term_wrong - C_term_wrong\n",
    "assert poly_big_wrong // T * T == poly_big_wrong\n",
    "## Error: A*B - C can not divide T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b512196f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A_random:  7\n",
      "B_random:  23\n",
      "C_random:  52\n",
      "Q_random:  64\n",
      "T_random:  53\n"
     ]
    }
   ],
   "source": [
    "## Linear PCP\n",
    "# random number generation, we choose 6 for education purpose\n",
    "r_random = GF(6)\n",
    "# check A, B, C\n",
    "A_random = A_term(r_random)\n",
    "B_random = B_term(r_random)\n",
    "C_random = C_term(r_random)\n",
    "Q_random = Q(r_random)\n",
    "T_random = T(r_random)\n",
    "\n",
    "print(\"A_random: \", A_random)\n",
    "print(\"B_random: \", B_random)\n",
    "print(\"C_random: \", C_random)\n",
    "print(\"Q_random: \", Q_random)\n",
    "print(\"T_random: \", T_random)\n",
    "\n",
    "assert A_random * B_random - C_random == Q_random * T_random\n",
    "# pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "cf662e4a",
   "metadata": {},
   "outputs": [
    {
     "ename": "AssertionError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m/var/folders/xz/sj6gs0f150n67pbx3mg32v3c0000gn/T/ipykernel_65295/743810614.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0mT_random\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mT\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr_random\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m \u001b[0;32massert\u001b[0m \u001b[0mA_wrong_random\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mB_wrong_random\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mC_wrong_random\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mQ_random\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mT_random\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     13\u001b[0m \u001b[0;31m# Error: wrong witness, A*B - C != Q * T\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAssertionError\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# Soundness: wrong witness\n",
    "## Linear PCP\n",
    "# random number generation, we choose 6 for education purpose\n",
    "r_random = GF(6)\n",
    "# check A, B, C\n",
    "A_wrong_random = A_term_wrong(r_random)\n",
    "B_wrong_random = B_term_wrong(r_random)\n",
    "C_wrong_random = C_term_wrong(r_random)\n",
    "Q_random = Q(r_random)\n",
    "T_random = T(r_random)\n",
    "\n",
    "print(\"A_wrong_random: \", A_wrong_random)\n",
    "print(\"B_wrong_random: \", B_wrong_random)\n",
    "print(\"C_wrong_random: \", C_wrong_random)\n",
    "print(\"Q_random: \", Q_random)\n",
    "print(\"T_random: \", T_random)\n",
    "\n",
    "assert A_wrong_random * B_wrong_random - C_wrong_random == Q_random * T_random\n",
    "# Error: wrong witness, A*B - C != Q * T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f27f3626",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GF(42, order=67)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_random * B_random - C_random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81bc6986",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
